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Section 1 

INTRODUCTION AND SUMMARY 

The major advantage foreseen in manufacturing products in space is the 
near -absence of gravity. It may eventually prove feasible, in this unique en- 
vironment, to produce products which are superior in quality to those made on 
earth. In manufacturing processes involving confined fluids which are heated, 
the low gravity environment of space should virtually eliminate natural fluid 
convection driven by gravity. However, there are driving forces other than 
gravity which covild possibly produce significant fluid circulation in confined 
systems. If a free liquid surface is present, surface tension gradients due to 
temperature variations can induce convection. Thermally induced expansions 
and compressions in a confined gas could also generate fluid circulation. Con- 
vection driven by mechanisms other than gravity has received very little at- 
tention to date; thus the current effort was initiated to determine the magnitude 
and effects of non-gravity convection on typical space manufacturing processes. 

The basic approach used to analyze this problem is to formxalate a math- 
ematical model of a simple yet representative system and obtain solutions for 
typical boundary conditions. The trends and magnitudes of non-gravity con- 
vection can then be estimated. The sample problem chosen is that of a single 
component compressible fluid in a low -gravity confined region which is heated 
along a wall. A one -dimensional flow situation is assumed and the governing 
partial differential equations are developed from the full equations for conserva- 
tion of mass, momentum and energy. Since these are highly nonlinear and 
strongly coupled equations, a numerical solution utilizing finite -differences is 
employed. An explicit finite -difference scheme is used to program the equations 
on a Univac 1108 digital computer. Sample problems are solved in which 
thermally induced wave motion and heat conduction are included. Results have 
been obtained which indicate that pressure and thermal expansion effects can 
induce significant fluid motion and increase heat transfer in low gravity. In 
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addition, the pressure waves are confirmed to be acoustical as discussed in 
the literature. From the study it is concluded that non- gravity -driven thermo - 
acoustic convective motion in a compressible fluid could have significant effects 
on space manufacturing processes involving heated fluids. A recommendation 
is given for developing a two-dimensional analytic model for further study of 
pressure and thermal expansion convection in low gravity. 

Section 2 of this report details the analytic formulation of the problem 
including the conservation equations, boundary conditions, dimensional analysis 
and assumptions. A discussion of the numerical solution method and finite- 
difference equations is given in Section 3. Appendixes A and B discuss the 
computer program which implements the numerical algorithms. Section 4 
presents numerical resxilts for two sample configurations — infinite parallel 
plates and concentric cylinders. Resxilts for a one -dimensional radial model 
of the Apollo 14 Heat Flow and Convection Demonstration are also given. Con- 
clusions are drawn from the res\ilts and recommendations for further analysis 
are given in Section 5. 
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Section 2 

ANALYTIC DEVELOPMENT 


2.1 THE PROBLEM 

In conjunction with NASA's Manufacturing in Space Program, several 
demonstration experiments have flown aboard Apollo spacecraft on their moon 
missions. The Heat Flow and Convection Demonstration Experiments (Ref. 1) 
were flown on the Apollo 14 mission to obtain data and information on heat 
transfer and fluid behavior in a low-gravity environment, A similar experi- 
ment package has recently flown on the Apollo 17 mission. The results of 

the Apollo 14 experiments strongly indicate that some mode of fluid convection 

-6 -4 

occurs even in a 10 to 10 g environment. The radial cell portion of these 
experiments exhibits a particularly interesting behavior. The "radial cell" 
consists of a cylinder of C02 gas which is heated at the center of the cylinder. 
Liquid crystals were used as temperature indicators. A theoretical heat 
transfer model of this unit was designed which includes conduction and radia- 
tion-but no fluid motion. A comparison of the theoretical model temperature 
predictions with the actual flight data is shown in Fig. 1 taken from Ref. 1. 
Although the magnitude of the numbers are in "reasonable" agreement, the 
shape of the flight data curve indicates that fluid convection was probably 
occurring. This curve rises faster at the early times than the pure conduc- 
tion curve and then levels off while the conduction curve continues to rise. 
Reference 1 provides more details on these experiments. 

The current analytic study of thermal convection was initiated to ex- 
plain the behavior of the Apollo 14 radial cell data and to develop the analytical 
capability necessary to assess the role of low-gravity convection in space 
manufacturing processes. 
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2.2 ANALYTICAL APPROACH 

The first step taken in this analysis was to postulate a possible driving 
mechanism for thermal convection of a completely confined compressible fluid 
in a low -gravity environment. The mechanism isolated for study here is the 
pressure and thermal expansion effects which result when a confined com- 
pressible fluid is heated. This problem, herein termed thermoacoustic con- 
vection, has received little attention in the literature. Trilling (Ref, 2), 

Knuds en (Ref, 3) and Luikov and Berkovsky (Ref. 4) used a linear perturbation 
analysis to investigate the wave motion induced in gases by boundary tempera- 
ture gradients. They found that sharp rises in boundary temperature can in- 
duce expansions which cause pressure waves to propagate through the fluid in 
much the same manner as pushing a piston through a gas filled pipe. Larkin 
(Ref. 5) investigated the thermal expansion effects in a confined gas in zero 
gravity. He solved the nonlinear conservation equations, with compressibility 
effects, using a finite -difference numerical scheme on a digital computer. His 
results indicate that heat transfer rates and pressure rises are significantly 
increased over predictions which neglect the thermally induced fluid motion. 

His calculations confirmed the acoustic nature of the velocity waves. More 
recently, Thursaisamy (Ref. 6) reached similar conclusions while investigating 
pressure behavior in spacecraft cryogenic tanks. 

The approach taken in this analysis is to formvilate a mathematical model 
of a simple yet representative fluid system and to study solutions for typical 
boundary conditions encountered in "space manufacturing." The fluid mechanics 
model begins with the conservation equations of motion, continuity, and energy 
in Euler ian coordinates. For a viscous, heat conducting Newtonian fluid, 
these equations are found in standard text books (Refs. 7 and 8) in vector form: 

Navier -Stokes Equation 

p ^ = F - VP - V X f jtx (X 7 X V)J + vf(X + 2p) v . vj (1) 
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Continuity Equation 


Se + p(v. V) =0 


( 2 ) 


Energy Equation 
Equation of State 

P = P{p, T) (4) 

These equations, with appropriate bovindary conditions, describe the 
flow and thermal behayior of the fluid. These are, of course, highly nonlinear 
and strongly coupled equations such that general solutions are not possible. 
Appropriate assixmptions and simplifications must be made if any solution to 
the convection problem is expected. Two models are used in this study to 
reduce the equations to a manageable form; however, the nonlinear behavior 
still renders analytic solutions impossible. A numerical s olution utilizing 
finite -differences is employed. 

2.3 MATHEMATICAL MODELS 

Two flow models are considered. Both incorporate the same basic 
assumptions, but differ in geometric aspects. Model 1 consists of two infinite 
parzillel plates which bound a compressible fluid. The flow sitxiation is thus 
one -dimensional and can be described in a Cartesian reference frame. Model 
2 consists of a radial segment of two concentric cylinders with a compressible 
gas between them. The outer cylinder then represents a fluid container and 
the inner cylinder is a heater. A one -dimensional radial coordinate system 
then describe the geometry. Both models utilize the following assumptions; 
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• Newtonian fluid obeying Stokes viscosity 
hypothesis {X = -2/3 fi) 

• Constant thermal properties k, C^, /i, y 

• No radiation or internal heat sources 

• No viscous dissipation of energy 

• Body forces are negligible (g/ge^<^l) 

• Ideal gas equation of state (P = pRT) 

Most previous investigators have invoked the classical Boussinesq 
approximation which neglects the effects of pressure on the density profile 
and allows density to vary with temperature only in the body force term. The 
density is thus constant in all other terms resulting in a quasi-incompressible 
approach. This assumption is not made here. The equations of compressible 
flow are used with variable density in all terms and related to pressure and 
temperature through an ideal gas equation of state. 

The primitive variable form of the equations are used as opposed to 
invoking transformations or combining equations to introduce fictitious variables 
The modern literature, (e.g., Ref. 9), indicates that this approach may yield 
results which are more physically correct and also may have advantages in 
numericcil stability and accuracy. The two mathematical models are now given. 

2,3.1 Infinite Parallel Plate Model 

The schematic below depicts the flow situation for Model 1. 

Isothermal Wall 

I 

X* = U — 1 


i 

t 1 


Compressible 

L' 

Fluid 

x' 

r 

i 



Heated Wall 





The goverhing equations (with the prescribed assumptions) can now be obtained 
from the general system (1) — (4), The one - dimensional Cartesian form is 
derived in most fluid mechanics texts (Refs, 7 and 8) and the derivation is not 
repeated here. Primes (’) are used to indicate that a quantity has physical 
dimensions; unprimed variables are dimensionless. The governing equations 
take the following form: 

Momentum 


d 

8t* 




3x' 


(p*u’u') 


-8 P* 
9x' 


+ 4/3j[l' 


u' 


3x' 


|2 


(5) 


Continuity 


8t» ^ ax' 


(p'u') = 0 


( 6 ) 






8T' 

8t« 


+ p'C^.u‘ 


8T' 

ax' 


-P' 


au« 

ax' 


+ k' 


8 ^T' 
8x'2 


State 


P' = p'R'T' 


(7) 


( 8 ) 


These eq\iations will be nondimens ionalized in Section 2,4 and expressed with 
dimensionless groups appearing as coefficients. 


The intial and boundary conditions for Model 1 are expressed mathe- 
matically as follows: 
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Velocity 


r 


u»{x',t«=0) = 0 

initially at rest 

u‘(x'=0, t‘) = u» (x‘=l!,t‘ ) = 0 

no slip at walls 

Temperature 


T‘(x‘,t‘=0) = Tq» 

initially isothermal 

T«(x'=0.t») = T^’(t') 

heated wall x* = 0 

T'(x'=L',t») - T^ 

isothermal wall x' = U 


(9) 


( 10 ) 


Pressure 


ap» 

Ox* 


t‘=:0 


= 0 


no body forces 


( 11 ) 


Density 


p»(xV,t»=0) 


Po' 

R’To* 


equation of state 


( 12 ) 


The thermal boundary condition at the heated wall is shown as a prescribed 
temperature history. This is done for simplicity of presentation. A pre- 
scribed heat flux boundary condition at x =0 can also be used in Model 1 as 
seen later in Section 4. 


Equations (5) through (12) formally define the mathematical Model 1 for 
analysis of thermoacoustic conyection in low grayity. 
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2.3.3 Concentric Cylinder Radial Model 


r 


The schematic below depicts the flow situation for Model 2. 



The governing equations for this configuration are again derived from the full 
system (1) — (4) (see Refs. 7 and 8). In terms of dimensional variables (') the 
equations take the following form. 


Momentum 


d 

8t* 


(p'u')+i, I7, (r'p’u'o') 


-8P« 

8r* 


4/3/x' 


8^u' 

+ ~ 

8u' u* 

8r'2 

r' 

8r* ”^T2 


(13) 


Continuity 



Jl J- 

r* 8r' 


(r'p'u*) = 0 


(14) 


Energy 


8T* 


8T* 


p'Cv'u-^ 


= -P' 


1 ^ 


8 

8r‘ 


(r'u*) 


+ k» 


8^T» 

8r»2 


+ 4 


8T» 

8r’ 


(15) 
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State 


P* = p'R'T' (16) 

The nondimensional form of these equations are actually used in the computa- 
tion as in Model 1. Dimensional analysis is discussed in Section 2.4. 

The initial and boundary conditions for Model 2 are expressed mathe- 
matically as follows; 

Velocity 


u'(r*,t*=0) = 0 initially at rest 

u'(r'=r-Q',t') = u'(r'=L',tf) = 0 no slip at walls 


(17) 


Temperature 


T» (r', t‘=0) = Tq* 
fT'(r*=L',t) = To* 


Ot« 

8r* 


= 0 

r'=L» 


T‘(r'=ro',t') = T *(t’) 


r*=rK* 


k‘ 


initially isothermal 

isothermal wall 
or 

adiabatic wall 


(18) 


prescribed temperature 
or 

prescribed heat flux 


(19) 


Pressure 


8P« 

8r* 


= 0 

t‘=0 


no body forces 


( 20 ) 
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Density 


P! 

p'(r',t'=0) -^p^t equation of state (21) 

Eqviations (13) through (21) formally define Model 2 for analysis of thermoacoustic 
convection in low gravity, 

2.4 DIMENSIONAL ANALYSIS 


Dimensional analysis is commonly used for a variety of scientific 
problems. If is basically a formal mathematical procedure which yields the 
dimensionless parameters associated with a particular problem. The pro- 
cedures for performing dimensional analysis vary considerably among authors. 
Ostrach (Ref. 10) and Heliums and Churchill (Ref, 11) present two approaches 
for dimensional analysis of natural convection problems. There are two major 
uncertainties in choosing reference values for natural convection probelms — 
the characteristic time and characteristic velocity of these systems are not 
always obvious. Since the present problem involves acoustic velocity waves, 
the isothermal sound velocity of the gas is used as the characteristic velocity. 
The characteristic time used is the time required for a wave to travel the 
length of the container. This same reference velocity is also obtained using the 
method of Ostrach (Ref. 10) by equating the inertial forces and pressure forces 
in the momentum equation. 


The following dimensionless variables are now introduced: 


X 




T» u* p P« 

To' Po' 



( 22 ) 


Model 1 is used here to illustrate the dimensionless form of the differ- 
ential eqxiations. Model 2 is very similar and is not shown here, but is used in 
dimensionless form in all computation. The expressions in Eq. (22) are 



f 


substituted into Eqs. (5) through (8) and the boundary conditions Eqs, 
(12), Algebraic manipulation then yields the following dimensionless 
and boundary conditions. 


Momentum (Dimens ionless) 


(pu) + -^ (P“U) = 




Continuity (Dimensionless) 


If + ^{pu) = 0 


Energy (Dimensionless) 


p|f + H = • pfx. + (sfVr) 


9^T 


8x' 


State (Dimensionless) 


P =pT 


Initial Conditions (Dimensionless) 

u(x, t=0) - 0 T(x, t=0) * 1 p(x, t=0) = 1 

Boundary Conditions (Dimensionles s ) 

u(xaO,t) = u(x5l,t) = 0 
T(x=0, t) = T^(t) 

T(x=l,t) = 1 


(9) through 
equations 


(23) 


(24) 


(25) 


(26) 


(27) 


(28) 
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The dimensionless groups which appear in the equations are: 


Prandtl Number Pr = 


M omentum Diff us ion 
Thermal Diffusion 


Reynolds Number Re = 


Pq Inertial Forces 


M' 


Viscous Forces 


(29) 


Peclet Number Pe * Re • Pr Convection 

Conduction 

C • 

Ratio of specific heats y = 


For the present problems of interest, y and Pr are of order 0(1) and Re 
is 0(1 0^). This dimensionless form of the equations is used in all numerical 
computation. 


This completes the formal development of the mathematical models used 
in this study of thermal convection in low gravity. The numerical solution 
algorithms are now presented. 
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Section 3 

NUMERICAL SOLUTION METHOD 
3.1 METHODS SURVEY 

The following brief review of the literature is not intended to be a complete 
survey; only the work more pertinent to the present problem is reviewed. Many 
excellent articles were examined during this study, but only those directly re- 
lated to the current problem of numerical computation of natural convection 
are presented here. 

Heliums and Churchill (Refs, 12 and 13) have applied finite difference 
calculations to natural convection for several types of problems. Explicit, 
time -dependent finite difference approximations were made to the conservation 
equations for mass, momentum and energy. Gravity was the only driving force 
considered. This work has pointed out that any method which is to be used on 
a computer should undergo a rigorous stability analysis and that the convection 
terms must be handled with much care. The Heliums -Churchill technique uses 
alternating forward and backward differences depending on the direction of the 
fluid flow. Computer storage and run time were found to be moderate. Agree- 
ment with previous results was excellent in some regions and good in others, 

Fromm (Ref. 14) presents some numerical results for the Benard problem 
of heating a fluid layer from below in a gravity field. The classical Boussinesq 
approximation is made and the governing equations are transformed to the 
vorticity -stream-function form. Central finite differences are used in the 
numerical solution algorithm. Cases examined include Rayleigh numbers from 

7 

the critical value to Ra = 10 . Correlations of heat transport with Rayleigh 
number are given and excellent comparisons are made with previous analysis 
and experiments. 
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Clark and Barakat (Ref, 15) have applied numerical computation to the 
problem of two-dimensional, laminar, transient, natural convection of a fluid 
in a rectangular container with a free stirface. Explicit finite difference 
methods were applied with success to the problem of a vapor -liquid interface. 
They concluded, however, that implicit methods may be preferred if only steady 
state results are needed. The theoretical results are compared to experimental 
measurements from the literature and indicate qualitative agreement. 

Wilkes and Churchill (Ref. 16) studied natural convection of a fluid con- 
tained in a long horizontal enclosure of rectangular cross section. Two- 
dimensional motion was assumed. The vorticity and energy equations were 
solved numerically by alternating direction finite difference methods. Com- 
pressibility effects were neglected, the density was allowed to vary only in 
the buoyancy term, and gravity was the only driving mechanism considered. 
Numericcil results for several cases (heating from the side) were compared 
to experiment with good agreement. Their results demonstrate that finite 
difference methods can adequately simvilate the thermal convection problem of 
heating from the side. 

Larkin (Ref. 17) gives a brief but important discussion of the numerical 
conservation principles which must be observed in solving the continuity 
equation. He presents an implicit technique which preserves an overall mass 
balance in a closed system. It is easy to implement due to the tridiagonal 
form of the coefficient matrix, 

Samuels and Churchill (Ref. 18) used finite difference methods to compute 
hydrodynamic instability due to convection in a horizontal rectangular region 
heated from below. Critical Rayleigh numbers were determined for a series of 
Prandtl numbers and length -to -height ratios. One of the major objectives of the 
work was to assess the usefulness of finite difference techniques for computation 
of natural convection. The governing equations were transformed such that a 
vorticity equation and an energy equation were treated and an implicit alternating 
direction finite difference method was used to solve the equations. Calculated 
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critical Rayleigh numbers were found to be in excellent agreement with other 
analytical results for Prandtl numbers greater than unity. For Prandtl numbers 
* less than unity, the calculated critical Rayleigh numbers exhibited a dependence 
on Prandtl number, a fact not predicted by linearized theory. If symmetrical 
initial conditions were imposed, a non -unique set of solutions were obtained. 
However, if an asymmetric initial condition was imposed, the calculation 
always converged to a single unique flow pattern, 

Aziz and Heliums (Ref, 19) present results of numerical calculation of 
three-dimensional laminar natural convection. The complete Navier -Stokes 
equations are transformed and expressed in terms of vorticity and a vector 
potential, A finite -difference method using alternating directions is used for 
the parabolic equations (vorticity and temperature) and a successive over- 
relaxation technique is applied to th® elliptic stream function equation. Com • 
parisons with previous works for two-dimensional problems are made and 
the conclusion is reached that the authors* method has important advantages 
in speed and accuracy. This is the most complete work that has been found on 
the three-dimensional natural convection problem, 

Torrance (Ref. 20) presents an excellent summary, review and compar- 
ison of finite -difference computations of natural convection. Five numerical 
methods, were compared for calculating twovdimensional transient natural 
convection in an enclosure. Both implicit and explicit procedures were con- 
sidered. Consideration was given to stability, accuracy and conservation of 
each method and it was concluded that no one method has all the features that 
are desirable. An explicit procedure developed by Torrance was shown to be 
conservative and stable without a restriction on the spatial mesh increment, 

A tentative conclusion reached by Torrance is that, (1) the Dufort-Frankel 
method will require less computer time if the mesh size restriction can be 
satisfied, (2) the Torrance fnethod should be used if the results are inter- 
preted in the sense of a large truncation error. The paper presents an ex- 
cellent comparison chart of the type of difference forms which various authors 
have used. 
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Plows (Ref. 21) presents numerical solutions for the laminar Be'nard 
problem. The Boussinesq-type equations are solved by a numerical scheme 
utilizing centered "leapfrog" differences for first order derivatives and a 
Dufort-Frankel pattern for second derivatives. An iterative scheme is used 
to determine Nusselt numbers and roll patterns for Rayleigh numbers between 
2000 and 22»000 for a range of Prandtl numbers. Nusselt number versus 
Rayleigh number calculations are presented and compare favorably with 
previous analytical and experimental studies, 

Schwab and DeWitt (Ref, 22) conducted a numerical investigation of free 
convection between two vertical coaxial cylinders. The ccMipled, nonlinear, 
partial differential equations were converted to finite -difference form and 
solved using an alternating direction implicit procedure, Restdts, in the form 
of steady state contour maps, are presented for several combinations of Prandtl 
and Grashof numbers. Among the conclusions reached is that a fully developed 
boundary layer flow was found to exist in the cavity for Rayleigh numbers 

3 

greater than 5x10. The variation of steady state Nusselt numbers with 
Prandtl and Rayleigh numbers and with geometric ratios was also briefly 
discussed. 

Cabelli and DeVahl Davis (Ref. 23) present a numerical study of the 
Benard cell problem for the case where buoyancy and surface tension are 
coupled. The conservation equations are expressed in the vorticity- stream 
function form and a surface tension boundary condition is imposed on the free 
surface. The vorticity and energy equation were solved with an alternating 
direction implicit finite -difference scheme and the stream-function equation 
was solved by over -relaxation. This appears to be the first work to solve the 
full conservation equations with surface tension effects coupled to buoyancy. 

The numerical results indicate that surface tension effects encourage the 
natural fluid motion in liquids. The reader is referred to this reference for 
details of their calciilations. 

Heinmiller (Ref. 24) has developed a mathematical model of thermal 
stratification and natural convection occurring in supercritical oxygen under 
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low -gravity environments. His model uses an explicit finite difference 
solution of the primitive equations for conservation of mass, momentum and 
energy. The Boussinesq assumption is not made and real gas properties are 
used throughout the computation. This work appears to be among the first to 
solve the fxill compressible form of the equations for a two-dimensional problem. 
The numerical model was used to successfully simulate a pressure collapse in 
the Apollo 12 oxygen storage system due to an acceleration change. 

Barton et al, (Ref, 25) also developed a numerical model for analysis 
of the Apollo spacecraft oxygen tank system, A numerical algorithm known as 
the General Elliptic Method (GEM) is developed for solution of the full conserva- 
tion equations in terms of the primitive variables. This method was also applied 
with success to the Apollo oxygen tank stratification problem. 

Thuraisamy (Ref. 6) presents a detailed one -dimensional model of 
natural convection in zero gravity. His aim was to analyze the flow of super- 
critical oxygen in the Apollo tanks. His model is a one -dimensional radial 
segment through two concentric cylinders - the outer cylinder being a tank 
wall and the inner cylinder being a heater. The model includes pressure and 
thermal expansion effects which are the subject of this report. Details of this 
work will be given in Section 4 where comparisons with present calculations 
are given. 

Larkin (Ref, 5) was the first to apply numerical computation to the 
governing equations including compressibility effects. He considered the one- 
dimensional flow of a perfect gas in a zero gravity confined container. Explicit 
finite differences were used on the momentum and energy equations and an 
implicit method was applied to the continuity equation. The results show that 
the pressure wave motion is acoustical in nature. Heating rates and pressure 
rise are greatly increased over those obtained with conduction alone, Larkin’s 
analysis will be detailed in Section 4 and comparisons will be made to the present 
work. 
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A study of other applicable numerical methods was made during this 
literature survey. Textbooks such as Refs, 26 and 27 are excellent for 
general problems but application to thermal convection is lacking. References 
9 and 28 through 30 contain a variety of finite difference techniques for applica- 
tion to general fluids problems, Cheng (Ref. 9) presents an excellent discussion 
of the conservation law approach to the numerical solution of the Navier -Stokes 
equations. Lax and Wendroff (Ref, 28) also discusses the conservation law 
form for solution of hyperbolic equations. Campbell (Ref, 29) presents a stability 
analysis for a finite difference solution to the Navier -Stokes equations, Brunson 
and Wellek (Ref, 30) present a mathematical discussion of the numerical stability 
of a Dufort-Frankel scheme for solving systems of parabolic equations, 

3.2 FINITE - DIFFERENCE EQUATIONS 

An explicit finite difference scheme will be applied to the present 
thermal convection problem. The unsteady form of the equations are used to 
study the transient behavior of the problem. The explicit approach with the 
unsteady equations allows a forward -marching -in -time algorithm to be employed. 
Forward time differences are used on the unsteady terms and space -centered 
differences are used on all space derivatives except the convection terms where 
a "flip-flop" forward -backward scheme is used. Although this scheme is only 
first-order, it should be siafficiently accurate for the qualitative result sought 
here. This was verified by comparing results to a second -order Dufort-Frankel 
algorithm (see Appendix B). The explicit formulation was chosen primarily 
because it is simple, easy to program and is non-iterative in nature. 

A node centered spatial mesh, shown in the schematic on the next page, 
is used to write the difference equations. In this grid, denotes space point i; 

X. = (i - l/2)Ax. The grid spacing Ax is constant. 
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Boundary Conditions 


x' = L'* 


i = k- 




Node Center 

Node Boundary 





X 


All quantities, (u, p, T,P), are defined at node centers, and differences are 
taken across a node using the "half -increment” quantities. This approach 
offers better conservation features than a scheme based on mixed node- 
center/node-boundary differences. A similar grid system was used by 
Heinmiller (Ref. 24), The half -increment quantities are defined at node 
boundaries by interpolation over node centered variables, i,e,, 


Pi+i = 2(Pi + Pi+i> 

^-1= + ^i_i) 


(32) 


Time derivatives are approximated by forward differences as 
follows; 

n+1 „n 
•^i ' -^i 

At 


8T 

8t 


( 33 ) 



r 


Where i denotes space point x., n denotes time t and n+1 denotes the new 

1 n 

time First order space derivatives except convection terms are 

approximated by central differences. 


8P 

8x- 



2Ax 


(34) 


Second order space derivatives are approximated by the usual centered 
form; 


i!u 

3x2 




n ~ n , .n 
u . , , - 2u . + u . , 
1+1 1 1-1 


{ Ax)2 


(35) 


The convection terms must be handled in a special manner in the explicit 
approach. If space-centered differences are used, the method becomes unstable 
regardless of the step sizes (Ref, 26). Heliums and Churchill (Ref. 12) and 
Ijarkin (Ref. 5) have used the following "flip-flop” method successfully. 



This forward -backward alternating direction scheme is used depending on the 
direction of the recirculating flow. Reference 12 shows that this method is 
conditionally stable. In the present work, the convective terms in the momentum 
and energy equations are differenced accordii^ to this scheme. 

The explicit differencing of the continuity equation requires special 
attention. As pointed out by Larkin (Ref, 17), strict mass conservation 
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laws must be obeyed numerically if meaningful results arc expected. However, 
the usual centered difference explicit method which would numerically conserve 
mass is unconditionally unstable. Larkin used an implicit method on this 
equation in order to bypass this restriction. In the present work, an explicit 
method similar to that of Heinmiller (Ref.24) was devised which is conditionally 
stable and conservative. The basis of the m,ethod is the use of the divergence 
form of the momentum and continuity equations. The momentum equation is 
solved first to yield the pu product at the new time point. This pu product at 
the new time is then central differenced in the continuity equation to yield the 
updated density p. The velocity is then computed from (pu)/p. 


Special forms of the difference equations are used at the nodes adjacent 
to a wall. This is necessary for accurate handling of boundary conditions. 
For example; 



node 1 






(37) 


Forms such as this are derived by computing node boundary points using 
linear interpolation and the differencing across the node using the known 
boundary conditions. 

A complete listing of the difference forms used in this work are shown in 
figs. 2, 3 and 4 for the momentum, continuity and energy equations respectively. 
In this figure, i denotes space node xj, n denotes time point tj^; (i=l and i=k 
are the nodes adjacent to a wall). The finite difference eqioations are shown for 
Model 1 only since those for Model 2 are very similar. 

A Dufort-Frankel numerical scheme was also developed during this 
study. Details are given in Appendix B. 
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Momentum Equation 
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Finite Difference Form 
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Fig. 2 — Finite Pifference Form of Momentum Equation 
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Continuity Equation 

■ft = ■ <p“> 


Finite Difference Form 
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Equation 




Finite Difference Form of Energy Equation 
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3.3 NUMERICAL STABILITY 

The explicit finite difference scheme just presented is a conditionally 
stable method. If certain criteria are met, the difference eqviations that are 
solved will approximate the solution of the differential equations. If the 
stability criteria are not met, unrealistic solutions or no solutions at all 
will result. The stability criterion of this method is a restriction on the size 
of the time increment At. The stability of the method is free of a spatial 
mesh (Ax) restriction and only accuracy dictates the size of Ax. 

The time step restriction for this method was determined by using the 
available literature and by experimenting numerically. This approach was 
necessary since a mathematical analysis of stability for nonlinear problems 
remains beyond the state of the art. For the present problem, in which 
pressure wave propagation contributes significantly, the time step should be 
determined by the time required for a pressure wave to propagate over a 
distance of one node length Ax, i.e.. 


At« < 



min 


(38) 


where c* is the velocity of sound. This is recognized as the familiar hyper- 
bolic limit known as the CFL condition (Ref. 26). Numerical experiment has 
shown this to be the governing stability criteria for the present problem. 

Since |u'|«c‘ and c' = R’T' for a perfect gas, we have 


A t> < 


Ax 

\/yR*T* 


min 


(39) 


In terms of the dimensionless variables (Section 2.4) we get the following 
stability requirement. 
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At<-^= (40) 

For moderate values of Ax, i.e., 20 mesh points. At is of the order of 0,01 

-4 

which corresponds to about 10 seconds of real time. Even for a machine 
of the capability of the Univac 1 1 08 this problem is not trivial in terms of run 
time. It may appear that an implicit procedure would relieve this time step 
restriction. Note, however, that it is still the physical phenomena of acoustic 
wave motion that restricts the time step size. An implicit method, using 
larger time steps, may bypass important physical behavior, 

3.4 SCALING PRINCIPLES 


A scaling procedure developed in Ref, 24 and used in Ref. 6 has been in- 
corporated into the present numerical algorithm to lessen the severe time step 
requirement. This procedure is now briefly outlined including its applicability 
to the present problems and limitations on the magnitude of scaling which is 
permissible. 


The basic problem arises because diffusion processes occur on a time 
scale much larger than that of acoustic wave motion. The purpose of scaling 
is to speed up certain of the physical processes occurring in the fluid without 
disturbing the thermodynamic state of the fluid itself. The following procedure 
is based on familar similarity laws of fluid mechanics (Ref. 8 ). The dimension- 
less groups which apply to the present problem are: 


Reynolds number Re = 


0 «L'u» 


u* C«‘ 

Prandtl number Pr = ^ P 

k» 

a'L‘ 

Nusselt number Nu = 

Mach number M = — 7 - 

c* 


( 41 ) 
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The objective is to increase the real time/computer time ratio by a factor 
of s such that each time step At that is used for computation will correspond to 
a real time s At, To achieve this objective, we define a flow situation, which is 
similar, but not perfectly equivalent, to our original problem, by the following 
tr ans f or mat ions : 


t’ = St' 

s 

k's = sk' (42) 

q'g = sq' 

u's = su' 

where the s subscript indicates a flow variable in a new time frame. We then 
have the thermal conductivity, viscosity, heat input and flow velocity increased 
by a factor of s. The new time frame ts in which we compute will now correspond 
to real time, st. The thermodynamic state of the fluid, temperature, pressure 
and density, remain unchanged in the new time system. We must now note 
that each of the dimensionless groups in Eq, (41), except Mach number, are 
the same in both time frames. We have increased the Mach number by a factor 
of s since the properties of the fluid which determine the sonic velocity have been 
preserved. Since the sonic velocity is uhchanged, the size of the permissible 
computation time step Ats is unchanged, but the real time At is increased. This 
satisfies our original objective. 

We can assess the limitations on this scaling procedure by considering 
the differential Eqs, (5) and (7), The scaling method essentially multiplies 
the unsteady terms and the diffusion terms by a factor s, or equivalently 
dividing the convection terms by s. The procedure is thus valid for problems 
in which the convection terms are "small” relative to the other terms in the 
equations. The size of the factor s is then a function of the "smallness" of these 
convection terms. For s = 1 we have an exact simulation of the problem and 
as s approaches infinity, we approach as a limit the pure conduction solution 
for a compressible fluid. 
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In the present analysis, calculations are made both for problems in 
which the convection terms dominate and for problems in which they are 
relatively small. The size of the scale factor s used in actual calculations 
is discussed in Section 4 where some results for typical problems are given. 

3.5 COMPUTER PROGRAM 

The finite difference method given in Section 3.2 was programmed in 
FORTRAN V language for a Univac 1108 Exec 8 digital computer. The 
equations for Model 1 and Model 2 were coded in the same program, termed 
TCI, with optional flags to control the program flow. A driver program which 
Cedis subroutines is used to allow complete flexibility of use. All input and 
output is performed in separate subroutines and each of the four basic equa- 
tions, momentum, continuity, energy and state, are coded in individual sub- 
routines. All data are input with punched cards and all output data are in 
printed page format. Provisions for running a problem in parts (timewise) 
is included. Details of the TCI program are given in Appendix A. 
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Section 4 

DISCUSSION OF RESULTS 
4. 1 SCOPE OF RESULTS 

This section summarizes the results which have been obtained to date 
for the one-dinnensional models of thermoacoustic convection. The objectives 
of the sample calculations were to 

• Qualitatively determine the importance of thermoacoustic 
effects as a heat transfer mechanism. 

• Assist in the analysis of the data from the Apollo 14 Heat 
Flow and Convection Demonstration. 

• Verify the numerical calculation method which may be 
used for analyzing other types of convection encountered 
in Space processing. 

These objectives are met by performing calculations, using the TCI program, 
for two representative problems. Case 1 consists of two infiiiite parallel plates 
with helium gas as the fluid. One wall is maintained isothermal and the other 
is subjected to an instantaneous temperature step. This case was chosen because; 
(1) solutions exist in the literature which can be used to verify the model and 
numerical method, and (2) this case should provide an upper limit on the 
severity of the thermal boundary conditions. Results for case 1 are compared 
and constrasted to those of Larkin (Ref. 5) and Thuraisamy (Ref. 6). The 
solution is shown to converge smoothly to an accurate steady steady given by 
an analytic expression. The importance of pressure convection is also 
established for this boundary condition. 

Case 2 consists of a cylindrical configuration with carbon dioxide as the 
fluid. The dimensions of the container and the fluid properties correspond to 
those of the Apollo 14 Heat Flow and Convection Demonstration (HFC). The 
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reader should refer to Ref. 1 for details of the HFC experiment. Three types 
of thermal boundary conditions are applied to case 2 — an instantaneous 
temperature step, a constant heat flux and a wall temperature>time history 
that corresponds approximately to the Apollo 14 HFC boundary condition. 

Steady state convergence is obtained for boundary condition 1 and compared 
to an analytic steady solution. Results for boundary condition 3 are compared 
to a pure conduction analysis of the HFC configuration. 

All of the results are shown in dimensionless units except HFC case 3 
solutions which are given in dimensional form in the International System of 
Units. 

4.2 PARALLEL PLATE MODEL 

The first case considered uses the Model 1 equations with helium gas 
as the fluid. The gas is initially at rest with uniform temperature. At time 
to, the X = 0 surface is instantaneously raised to twice the initial temperature. 
The problem configuration and fluid properties are illustrated in Fig. 5. 

Larkin (Ref. 5) and Thuraisamy (Ref. 6) have also applied finite -difference 
calculations to this problem and comparisons to their solutions are made 
whenever possible. 

The steady state solution to this problem can be obtained analytically. 

At steady state the fluid motion will vanish; thus we have 

= 0 (43) 

where the subscript indicates a steady state value. The temperature profile 
should then approach the solution of the conduction equation with the prescribed 
boundary values. This solutibn is 

T = 2 - 3C (44) 

ss 
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Property 


Value 



15,3 cm 

1.9 X 10”^ gm/cm^ 

1.01 X 10^ dyne/cm^ 

273 °K 

1.875 X 10”^ gm/cm-sec 
0.783 cal/gm-°K 

3.44 X 10"^ cal/cm-sec-^K 

1.66 

7.58 X 10^ cm/sec 
2 To 
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From the momentum equation we then require 

= 0 or P _ = constant 
ss 

The value of this steady pressure can be found from the ideal gas equation of 
state by using the mass conservation condition that 


ax 



pdx = 1 


We can integrate the p = P/T profile to obtain 


ss 


1 

In 2 


::: i-44 


and 


ss 


1.44 

2 - X 


(45) 


(46) 


A scale Victor a = 1 was used to calculate the transient solution to this 
problem to achieve an exact numerical simulation of the problem. The unsteady 
solution profiles are illustrated in Figs. 6-8. The sudden change in the wall 
temperature at x =0 induces an expansitm of the gas due to its compressibility. 
This expansion creates fluid motion in the +x direction and the waves propagate 
to the X = JL surface and then reflect. The subsequent velocity field is then 
influenced by the reflections of this disturbance from each of the plates. 

Figure 6 shows the calculated dimensionless velocity profile as a function 
of time for the first 400 time steps (At = 0.025). From this figure the oscillatory 
nature of the wave motion is seen. The period of the calculated waves is 1.55 
units of dimensionless time. The acoustic wave period in this system of 
dimensionless units is 2,/y/y and for helium is 1.55 units. The calculated waves 
are thus acoustic. This was also confirmed by Larkin’s (Ref. 5) calculations. 
The maximum amplitude of these first few waves is approximately 0.02 units. 
Larkin's calculation yields an amplitude of approximately 0.2 units with a period 
of 1.55 units. The present calculation thus yields an order -of- magnitude lower 
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amplitude than previously predicted. Thuraisamy (Ref. 6) used Larkin's 
numerical method to solve this problem and also calculated 0.2 for the amplitude. 

Private communication with both Larkin and Thuraisamy and with S. W. 
Churchill of the University of Pennsylvania has resulted in the conclusion that 
the differences in the numerical method of Larkin and the present method are 
probably responsible for the note^ discrepancies. The present calculation of 

i 

0.02 units appears to be more correct for the following reasons. 

The quantify u /^RT^ is a Mach number and would be expected to be 
"small" (relative to Mach 1) for natural convection flows. A Mach number of 
0 02 appears to be more physically realistic for the present problems than 0.2. 

This is further confirmed by analyses bn the similar problem of Sondhauss 
tubes by Feldman (Ref. 32) which indicated that the Mach number was always 
less than .06. To checkout this discrepancy’ further, a second order finite- 
difference method -- The Dufort Frankel scheme ^ was developed at Lockheed 
independently of the present method. As discussed in Appendix B, this Dufort 
Frankel method also produced a calculation of approximately 0.02 for the wave 
amplitude. For these reasons the results of the present analysis shown in Fig. 6 
are considered to be the- true solution for the velocity profile . 

The calculated temperature profile at t = 1000 units is shown in Fig. 7. 

The present calculation is in excellent agreement with those of Larkin and 
Thuraisamy. The solution at t = 1000. is seen to-be much closer tO the steady 
state than a conduction solution which neglects thermally induced fluid motion. 

This figure indicates that thermoacoustic convection is an effective transfer 
mechanism and greatly inhances the rate of heat transfer. It is interesting to 
note that the temperature profiles calculated by the present method are in 
excellent agreement with those of Larkin while the velocity amplitude calcu- 
lations are quite different. 

Figure 8 shows the calculated density profile and mean pressure distribution 
for this problem. The density distribution at t = 1000 has almost achieved its 
steady state profile. The spatially averaged pressure versus time is compared 
to Larkin's result. The agreement is excellent with only a slight discrepancy 
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Fig. 8 - Density and Pressure Profiles for Infinite Plate 
Problem (Helium, T =2) 







at the early times. Note the large increase in pressure over that predicted 
by a Rayleigh-type solution. The Rayleigh solution is obtained by dropping 
all convective terms in the energy equation and computing mean pressure 
from the integral of the temperature profile with density constant, i.e. no 
convection. The strong influence of thermoacoustic convection on the rate of 
pressure rise is clearly seen. It is interesting to note that t = 1000 units 
corresponds to approximately 0,2 second of real time. The 30% rise in 
pressure in this short time may appear physically unrealistic, however, we 
should recall the severity of the thermal boundary condition (T^ = 2T^). This 
result thus represents an upper limit on the rate of pressure rise due to 
the r moac ou stic ef f ect s , 

The calculations thus far have established thp>t thermoacoustic convection 
can greatly enhance the rates of heat transfer and pressure rise. Reasonable 
agreement with previous solutions for the unsteady flow profiles has also been 
illustrated. To verify the steady state convergence of the method, the helium 
problem was solved with a scale factor s = lOO to achieve a "large" number of 
time step with a reasonable amount of computer time. With this scale factor, 
minute details of the unsteady profiles will be missed but the steady state 
should not be effected. Figure 9 illustrates the convergence of the computation. 
A plot of Nusselt number versus time shows that a steady state (Nu = 1) is 
approached. The definition of the unsteady Nusselt number used here is 

8 T 

Nu = 8 X w (47 ) 



where 8 Tl is the dimensionless temperature gradient at the heated wall 
w 

(x = 0) and is the steady state temperature difference T^ T^. The 

figure shows that the temperature profile near the x = 0 wall approaches its 
steady state in approximately 1000 units of time. 

The mean pressure profile also smoothly approaches the true steady 
value of 1.44. Approximately 200,000 time steps or 50 seconds of real time 
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was required to reach this constant steady state pressure. The spatial pressure 
distribution however remained constant after 80,000 steps, tljus we have that 

8 P ► 0. 

8 X 

-8 

The velocity profile calculated after 80,000 steps showed a mean of ~ 10 and 

-12 

after 200,000 steps was of the order of 10 , i. e. the wave motion is damped 

out as expected. Figure 9 also shows the unsteady temperature profiles for 
three locations in the flow. These also remain unchanged after about 200,000 
calculation steps or 50 seconds. 

The preceding calculations were carried out using 20 mesh points and a 
time of 0.025 units. The numerical solution of this problem for the unsteady 
profiles (Figs, 6, 7 and 8) required 3 minutes of computer CPU time on 
the Univac 1108 computer using a scale of s = 1 for exact simulation. 

The steady state convergence (200,000 time steps) using s = 100 required 15 
minutes of CPU time. These computer time/real time ratios are not small 
but are well within the practical limits of the current state of the art. More 
sophisticated time-scaling methods could improve the computational economics 
even more. 

*. 

Several interesting results were found concerning assumptions which 
are classically made in the numerical computation of natural convection. These 
were explored to determine the importance of viscous dissipation and pressure 
convection on the solutions of the energy equation. The viscous dissipation 
term in the energy equation is usually neglected in numerical solutions of 
natural convection. An order of magnitude analysis was performed and it 
appeared justifiable to neglect this term. To verify this assumption, the 
TCI program was mo<|ified!to include the dissipation term 



The infinite parallel plate problem was solved With <j> included and compared to 
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the solution neglecting dissipation. Typical results at t = 1000 are shown below. 
This clearly verifies the relative unimportance of Viscous dissipation even 
for a severe thermal boundary condition T„, s 2. 

Vv 


Quantity 

Without 4 

With <f> 

T (X s .025) 

1.9603 

1.9611 

T (x = .975) 

1.02175 

1.01302 

p ( x= .025) 

0.6749 

0.6751 

p(x=.975) 

1.305 

1.306 

P 

1.3266 

1.3261 

avg 



Nu 

i 

1.286 

1.254 


Classical formulations of the thermal convection problem have utilized 
the Boussinesq approach which treats quasi -incompressible fluids as discussed 

in Section 2. Most studies which have considered compressible flows have not 
given full consideration to the effects of spatial pressure gradients on the flow 
work in the conservation of energy equation. A simple, yet infornaative study 
was made to determine the consequences of such an approximation on the un- 
steady flow profiles. 

Consider the energy equation Eq. (7) written in terms of the constant 

pressure specific heat C . The flow work term for one dimensional geometry 

P 

is then 


9P + u 3P 

at 3x . (49) 

The present analysis to this point has included the full flow work term Eq. (49). 
For the purpose of evaluating the importance of the pressure convection term 

a P 

u , the following modification was made. 
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d P 

• The u Is assumed small and is neglected. 

8 P 

• is calculated from a spatial mean pressure and 
varies only vsrith time. 

This approach follows a classical method of analysis for natural convection. 

This case was solved using the helium parallel plate problem to compare the 
solutions with and without pressure convection. 

Figure 10 shows profiles of pressure and temperature versus time with 
and without the u 8 P/8x pressure convection term. It is evident from these 
curves that the term cannot be neglected for the sample problem studied. 

There is considerable deviation in the unsteady profiles for both temperature 
and pressure. It is interesting to note that the temperature is higher without 
the pressure term than with it. Neglect of 8P/8x implies that the pressure is 
constant across the region between the two plates. Since there is a thermal 
gradient across the fluid field, this implies that the pressure equalizes 
instantaneously; i. e., the pressure waves propagate at an infinite speed rather 
than at the true speed of sound. This is shown by Fig. 10 since this figure 
indicates that heat from the wall is being transferred to the fluid faster due 
to the infinite wave speed. Note that both the pressure and temperature profiles 
achieve their steady state much more rapidly if the finite wave speed is neglected. 

For this sample problem we must conclude that neglect of pressure con- 
vection is not justified. Although no general conclusion can be made from this 
one case, we are alerted to the fact that the finite wave speed may be important 
in studying convection of fluids which are rapidly heated. A more detailed 
parametric study, with heating rate as the parameter, could be performed to 
determine the ranges of heating rate where this effect is important. 

4. 3 RADIAL MODEL 

Computations are presented for a second configuration which uses the 
Model 2 equations with carbon dioxide as the fluid. This model represents a 
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one -dimensional radial segment of the Apollo 14 Heat Flow and Convection 
Demonstration (Ref. 1) radial cell experiment. The HFC radial cell, the one- 
dimensional mathematical model, and a list of fluid properties used, is given 
in Fig. 11. Three types of thermal boundary conditions are applied to this 
configuration: 

• Instantaneous wall temperature step 

T ' = 2T ' , 
w o ’ 

• Constant heat flux 
q^' s constant 

• Wall temperature-time history prescribed 
T^' = f ' (f ) 

These are boundary conditions on the r = r_ surface, i. e., the inner cylinder 
which is a heater. The outer wall is held isothermal except in the constant 
heat flux case where an adiabatic outer wall is assumed. The steady state 
solution for this model using the T^’ = T^' boundary condition can be easily 

established analytically but not for the other boundary values. As with Model 1, 
the "one -dimensional" motion will vanish at steady state, thus 

u„ = 0 (50) 


where the subscript indicates a steady state value. The temperature profile 
will then approach the solution of the conduction equation with the prescribed 
boundary conditions. This solution is 



In(r^r) 
In r 

o 


( 51 ) 
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From the momentum equation, we then require 


8P = 0 
8 r 

ss 


or P = constant 
ss 


The value of this steady pressure can be found from the ideal gas equation of 
state by using the mass conservation condition that 



We can integrate the p - Tj, profile to obtain an infinite series for the steady 
steady state pressure. Taking r^ as in Fig. 11, we get 

Pss~l.l6 (52) 


and 


P 


ss 


1.16 In r 

o 

In (r r ) 

' o ' 


(53) 


A scale factor 9 = 1 was used for the case T = 2T to examine the 

w o 

exact nature of the wave motion and profiles for the cylindrical model. The 
unsteady flow profiles for this case are given in Figs. 12 through 14. The 
oscillatory velocity profile shown in Fig. 12 is similar to that produced for the 
parallel plate model (Fig. 6), The period of the acoustic waves is 1.62 with 
a maximum amplitude of 0.03 units. The wave amplitude is slightly higher 
than the plate model and the waves themselves appear more "irregular" than 
those in Fig. 6. The amplitude is still of the order of 0.02 and not 0.2 as dis- 
cussed in Section 4.2. 
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Fig. 12 — Velocity vs Time at r = .571 for Radial Model (T 


r 


Figure 13 gives temperature, pressure and density profiles for the 

I I 

= 2T^ case. The calculated temperature profile versus radial distance 
at t = 350 units again indicates the large increase in heat transfer compared 
to a conduction only model. The thermoa :oustic effect is also an important 
heat transfer mechanism for the radial configuration. The solution profile is 
seen to be approaching a true steady state condition which it achieves near 
t s 1000 units. The spatial mean pressure rise versus time is again con- 
siderably above that of a conduction (Rayleigh-type) solution but is not as 
drastic a rise as given by the plate model (Fig. 8), The density profile is 
shown as a function of radial distance at t s 1000 and indicates the smoothness 
of the steady state profile. These calculations are in excellent agreement with 
the analytical steady state given by Eqs. (50) thr ough (53. ). 

Figure 14 shows the spatial velocity distribution for t = 10 and t = 500. 
The smoothness of the profiles versus r indicates that the numerical method 
is very satisfactory for calculating the relatively slow flow of natural con- 
vection. This case was also run for a scale s = 100 and profiles were calcu- 
lated to t = 100,000. The damping out of the velocity profile gives good indica- 
tion of the convergence of the numerical method. The velocity continued to 

-14 

damp out with tinie to a value about 10 , The temperature, pressure and 

density profiles at t s 100,000, or about 2 seconds, agre4 with the analytical 
steady state to less than 0.1% deviation. 

The computation time for this case was 3 minutes on the Univac 1108 
Exec 8 using 20 mesh points and a time step At = .025. 

The next case for the radial configuration consists of applying a constant 
heat flux boundary condition to the r = r^ surface. A dimensionless heat input 

a ' L' 

q = gr^. = 20 (54) 

was used as the boundary value. When applied to the surface area of the inner 
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Radial Distance, r 



Radial Distance, r 


Fig. 14 - Velocity vs Radial Distance r for Two Time 
Points (Radial Problem T = 2) 
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cylinder, this corresponds to a heat input of 5,7 watts which is the heater 
power used to drive the HFC radial cell experiment. The r = 1 surfece was 
considered adiabatic for this case. 

A scale fector s = l was used to calculate the unsteady flow profiles to 
t = 2000 units. The velocity waves and pressure rise for this constant heat 
flux case were much less severe than the T ' s 2T ' boundary condition. 

However, the unsteady temperature profiles for fluid points near the heated 
surface exhibit an interesting behavior as illustrated in Fig. 15. The con- 
vection solution is compared to a conduction only solution at r = r^ (at the 
heated surface) and r = 0.093 (near the heated surface). At the early times, 
the fluid at and near the heated wall gets hotter than by pure conduction in- 
dicating an effective fluid convective mechanism. At about t = 1200 the conduction 
solution has "caught-up" and the local temperature calculated by conduction 
alone at these r -locations continue to rise with time. The convection solution, 
however, tends to steady out with time at these same r-locations since the 
convection mechanism is now transfering the incoming heat to further portions 
of the fluid. This same behavior was exhibited at locations further from the 
wall at later times. Although no direct comparison with flight data is possible 
for this case, the behavior is qualitatively the same as that obtained in the HFC 
flight experiment, Fig. 1, Computation time was six minutes. 

A third type of boundary condition was applied to the radial configuration. 
This boundary condition consists of a profile of the r = r^ wall temperature 
as a function of time. The computer program was modified to perform a table 
look-up and interpolation at each time step to obtain the psuedo-constant wall 
temperature. The profile that was used was obtained from a detailed conduction 
network analysis of the actual HFC radial cell experiment. An inhouse Lockheed 
thermal analysis program (Ref. 31) was used to calculate temperatures at 100 
points in the HFC radial cell including the inner cylinder heater post. This 
heater post temperature profile was then used as an input bo\indary condition 
to the radial model convection program. The outer wall temperature history 
as calculated by the conduction network showed this surface to be essentially 
isothermal. These boundary conditions are shown graphically in the top curve 
of Fig. 16. 
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Fig. 16 - Heater Temperature and Velocity Profiles for Apollo 14 
HFC Radial Simulation (T^ vs t Boundary Condition) 
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A complete analysis of the HFC requires that 600 seconds of real time 
be simulated. For this reason a scale factor s = 100 was used for the calcu- 
lation of this model. This time scaling should yield a quiet accurate simulation 
due to the relatively slow heater temperature rise. Figure 16 shows the cal- 
culated spatial velocity distribution at t = 10 seconds and t s 600 seconds. At 
the early times, flow velocities of the order of 0.5 cm/sec are calculated. As 
the simulation time increases these velocities are damped out to values the 
order of 0.1 cm/sec at t = 600 seconds. It is interesting to note that the peak 
velocities occur at r'/L’ =0.4 to 0.5, i.e. , near the center of the radial 
segment. This same behavior was exhibited at the intermediate time points 
also. The velocity wave motion for this case was much more highly damped 
than those obtained using the severe ^ *= 2T^’ condition.. 

Figure 17 gives profiles of the calculated temperature versus time at 
radial locations r* = . 69 cm and r* = 1.65 cm. It should be noted at the outset 
that no flight data are shown on this figure for the following reason. The 
Apollo 14 HFC radial cell experiment uses liquid crystal temperature indicators 
which must be placed in view of a Data Acquisition Camera. This prohibited 
the acquisition of data in the CO^ gss at radial locations from the heater post. 
The one-dimensional radial convection model discussed here is used to obtain 
qualitative information related to the actual HFC demonstration. A two- 
dimensional analysis, discussed in Section 5, is needed in order to make com- 
parisons to actual flight data. 

However, Fig. 17 does show a very feimiliar behavior. The calculated 
temperatures are compared to a one -dimensional (radial) conduction only 
solution using the same numerical method for the energy equation as used for 
the convection analysis. Also shown for comparison are the results of the 
detailed two-dimensional (axisymmetric) conduction network analysis using 
the Lockheed conduction program (Ref. 31). At both the radial locations shown 
the convection solution profile rises higher initially than pure conduction and 
then ’’flattens" out. At t = 600 seconds, the temperatures are lower than those 
predicted by pure conduction. 
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This is the type of behavior exhibited by the flight data comparison 
given in Fig. 1 and discussed in Ref. 1. The magnitudes of the temperature 

differences shown in Fig. 17 are not as large as the HFC flight data analysis 
of Ref. 1 predicts. The largest difference in Fig. 17 is the order of 5% between 
convection and pure conduction. It is interesting to note the large deviation 
between the one -dimensional and two-dimensional conduction solutions. This 
suggests that multi -dimensional convection may also be important in determ- 
ining the proper magnitude of the temperatures for this configuration. It is 
also emphasized that gravity convection was not included in the solutions shown 
in Fig. 17. The effects of gravity-driven convection would produce a larger 
deviation between the conduction and convection solutions of Fig. 17, Gravity 
convection plus thermoacoustic effects could possibly explain the behavior of 
the Apollo 14 HFC data. 

The results of this model, displayed in Fig. 17, leads us to the following 
conclusions- 


• Thermoacoustic convection effects can qualitatively 
cause the type of behavior seen in the Apollo 14 
HFC radial cell experiment 

• The magnitude of the temperature differences between 
convection prof iles and conduction- only profiles are 
not as large as indicated by the flight data analysis of 
Fig. 1 

• Multi- dimensional convection effects may be important 
for the HFC configuration and 

• The coupling of gravity-convection with thermoacoustic 
effects could possibly explain the behavior of the HFC 
flight data. 
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Section 5 

CONCLUSIONS AND RECOMMENDATIONS 

5.1 CONCLUSIONS 

Several conclusions are evident from this numerical study of thermoacoustic 
convection in low gravity. These are summarized below in two categories — 
thermoacoustic effects and the numerical method, 

5.1.1 Thermoacoustic Effects 

The sample calculations presented in Section 4 have shown that thermo- 
acoustic convection can be an important heat transfer mechanism. Specifically, 
it was shown that: 


• The thermally induced wave motion is acoustical, 

• Thermoacoustic convection can greatly enhance the rate of 
heat transfer. 

• The mean pressure rise in a confined fluid ia more rapid 
due to thermoacoustic effects. 

• Pressure convection effects on the conservation of energy 
can strongly influence the transient fluid behavior. 

• The magnitude of the effects of thermoacoustic convection 
is a strong function of the severity of the thermal boundary 
conditions. 

• Thermoacoustic convection can cause flow phenomena 
similar to that observed in the Apollo 14 HFC flight 
demonstration. 

• Quantitative comparison of theory to flight data will require 
a multi- dimensional analysis of thermoacoustic convection. 


The general conclusion is summarized as follows j For low-gravity space 
processii^ situations involving confined compressible fluids which are heated, 
the thermally induced fluid motion should not be neglected in performing 
analytical design studies. 
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5.1,2 Numerical Method 


The numerical calculation method us ed in this study has proved to be 
very satisfactory for computation of natural convection flows. It has been 
shown that the method; 


• Is conditionally stable with the hyperbolic limit restricting 
the time step size. 

• Is sufficiently accurate for space processing applications. 

• Converges to an accurate steady state for the problems 
studied. 

• Is computationally economic by state-of-the-art standards. 

• Is readily adaptable to time scaling laws to reduce the 
computer time/real time ratio. 

• Should be applicable to multidimensional natural convection 
computation. 

• Can be readily adapted to the computation of other types of 
convection such as gravity driven or surface tension driven 
flows. 

The numerical computation technique derived in this study has a wide variety 
of applications in the analysis of convection in space manufacturing processes. 


5.2 RECOMMENDATIONS 


Recommendations resulting from this study are offered for further in- 
vestigation of thermoacoustic convection for space processing applications. It Is 
recommended that the future effort consist of two tasks conducted in parallel. 


Task 1. Develop a two-dimensional/axisymmetric computer 
program for analysis of thermal convection of com- 
pressible fluids. 

Task 2. Conduct a simple ground-based laboratory experi- 
ment to investigate thermoacoustic convection in 
confined fluids. 


The two-dimensional program shoxild include, as a minimum, gravity driven 
convection and thermoacoustic convection. The direction of the gravity 
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vector with respect to the heated surface should be arbitrary. Both rectangular 
and cylindrical geometries should be included as options. The proven numerical 
method used for this one -dimensional study can be extended for use inthe two- 
dimensional program. This capability will allow the coupling effects f f gravity 
and thermoacoustic waves to be analyzed. The program will also provide a 
base for further extension to analyze convection in space manvif actur ing 
processes such as crystal growth, melting and/or casting of metals and 
chemical separation processes. 

The experimental program should be kept relatively simple while re- 
taining enough sophistication to obtain quantitative data. The experimental 
apparatus should primarily consist of a container of gas, a heater capable of 
rapidly raising the gas temperature, a system of sensitive pressure trans- 
ducers, thermocouples, and recording instrumentation. The primary purpose 
should be to detect the pressure waves in order to study their amplitude, to 
measure gas temperatures for comparison with theory and to investigate the 
coupling of gravity to thermoacoustic effects. 

This parallel effort will pr ovide additional insight into the thermoacoustic 
phenomena, provide experimental data for comparison with theory, and will 
provide a basic convection computer program for analysis of future space 
manufacturing systems. 
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Appendix A 


A.l PROGRAM LOGIC 

The TCI program consists of three main blocks: 

Input Reads and prints all input data and 

initializes parameters 

Calculation Solves momentum, continuity, energy 

and state equations 

Output Prints convection solutions. 

A block diagram of the program is shown in Fig. A-1. 

The program is designed in subroutine form with a DRIVER program to 
control the subprogram calls. Figure A-2 is a flow chart of the TCI DRIVER 
routine. The coding is entirely in standard third generation FORTRAN. It has 
been run on the Univac 1108 Exec 8 system, bu-t can operate on any system 
having a FORTRAN IV or V compiler and at least lOK of core memory. No 
external storage devices are needed. 

Problems using Model 1 or Model 2 (as described in Section 2) can be 
run with the TCI program. The physical input data can be either in metric or 
English units. The program calculations are all performed in terms of dimen- 
sionless quantities (see Section 2) and the output is also in terms of dimension- 
less variables. All input is via cards which are punched according to the input 
guide which follows . 

A. 2 INPUT GUIDJ 

All input to the TCI program is via punched cards. No tapes or drums 
are used. The input primarily consists of the following card groups. 
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Read and print all input data. 


Solve momentum equation 


Solve continuity equation 


Solve energy equation 


Calculate pressure (ideal gas). 


Print out convection parameters. 


Loop to march forward in time. 


Go run another case. 


Fig. A- 2 - Flow Chart of DRIVER Program for TCI 
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Card 1 - Title card 

2 - Program control flags 

3 - Physical parameters and initial conditions 

4 - Thermal properties 

5 - Bovuidary conditio|is 

In addition, there are other optional cards as folldws: 


Card groups 6-7 
Card groups 8-11 
Card group 12 


Variable boundary condition tables 

Restart information 

Normal run termination card. 


Details of the input data, formats and options are now given. 


CARD 1 


Format (80A1) 


Variable Description 

TITLE Problem description card. Any FORTRAN characters 

can be used to identify the problem being run. 

CARD 2 Format (1415) 


Variable 


Description 


NN 

NPRNT 
1S T AR T 


ICOORD 


Number of nodes to be used In finite-difference grid 
(50 maximum, nominally 20). 

Print control flag. Program will print out every 
NPRNT time steps. 

Restart flag . If ISTART = 1, program will run from 
the standard initial conditions. If ISTART = 2, profiles 
of temperature, density and velocity at a given time are 
read -in and used as initial conditions. This option 
allows restarting a probleni if more time steps are 
needed or if program tetminates due to "max-time." 

Coordinate control flag . If ICOORD = 1 , the inifinite 
flat plate model equations are used. If ICOORD = 2, 
the radial coordinate — concentric cylinder model 
equations are used. 
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CARD 2 (Cont’d) 

Variable 

lUNITS 

ITYPE 

CARD 3 
Variable 
XL 
TZ 

RHOZ 

PZ 

DTIME 

TIMEMX 

RZ 

SCALE 


Format (1415) 

Description , 

Units control flag. If lUNITS = 1, metric units are 
assumed to be input. If lUNITS ^ Z, Ehglish units are 
assumed. The appropriate input units are shown with 
Cards 3, 4 and 7. 

Flag to control the program flow. ITYPE = 1 rxms the 
program in standard mode; the momentum, mass and 
energy equations are solved in full. ITYPE = 2 signals 
a conduction-only run. The energy equation with no 
convective terms is solved. This option can be used to 
compare the effects of convection on the temperature 
profile. 


Format (8E10.0) 


Description 

Distance L between the two plates (if ICOORD = 1) or 
radius L of outer cylinder (if ICOORD = 2), (cm or ft) 

Initial isothermal absolute temperature T of fluid at 
rest, (°K or *^R). ® 

Density of fluid at temperature T^ and pressure P^, 

(gm/cm^ or Ib/ft^). 

Initial pressure P^ (dyne/cm^ or Ib/ft^). Either P^ 

or can be input and the program will calculate the 

other from the ideal gas state equation. 

Dimensionless time step At. (At = At^V^R^^^/ij ) used 

to start the calculation. The program monitors the 
stability criterion and adjusts At to maintain stability. 

Dimensionless time at which the calculation terminates. 

(t = t' Vr^/l'). 

' max max * o' 

Radius of inner cylinder (heater) if 'ICOORD = 2. RZ 
should be zero if ICOORD = 1 . 

Scale factor, s, (dimensionless) used to scale problem 
•'time" to ease the small time step restriction (see 
Section 3.4). 


A-5 



r 


CARD 4 

Variable 

CONDZ 

VISCZ 

CVZ 

GAMMAZ 

Z 

CARD 5 

Variable 

IBOUND 


TW 

QIN 

ITOP 


Format (5E10.0) 


Description 


Thermal conductivity, k, of fluid, (dyne/sec°K or 
Btu/ft-sec-°R). 

Dynamic viscosity, fi, of fluid, (gm/cm-sec or Ib/ft/sec) 

Specific heat at constant volume, C , of fluid (dyne -cm/ 
gm-°K or Btu/lb-®R). ^ 

Ratio of specific heats, Y (dimensionless). 

Compressibility coefficient in equation of state (dimen- 
sionless), The present program can use any 0 < Z < 1 
in computing pressure from P = pZRT 


Format (15, 2E10.0, 15) 


Description 


Control flag for boundary condition on x « 0 (or r = rQ) 
surface. IBOUND s +1 a specified temperature T^ is 
used. IBOUND = +2 a specified heat flux q is used. 
Positive values of IBOUND will cause a constant Tw or 
q;^ to be used. A negative IBOUND allows tables of T^ 
or qw versus time to be read-in and the boundary 
condition will be variable as a function of time (see 
Cards 6,7), 

Dimensionless wall temperature to be applied to x = 0 
(or r = ro) surface. If IBOUND = +1 this TW will be 
used as a constant bovmdary condition; If IBOUND = -1 
TW is computed by linear interpolation from the input 
table (Card 7). 

Dimensionless heat flux (q^L^ /k^T^) to be applied to 

x = 0 (or r = r^) surface. If IBOUND = +2, this QIN 

will be used as a constant boundary condition; if 
IBOUND = -2, QIN is computed by linear interpolation 
from the input table (Card 7). 

Boundary condition control flag for x = L surface. 

ITOP = 1, an isothermal surface (T ) at x = L will be 

o 

used. If ITOP = 2 an adiabatic surface (q = 0) at x = L 
will be used . 
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CARD 6 (Optional) Format (15) 

Variable Description 

NPTI Number of pairs to be input in T versus time or 

versus time tables. Input onl/ if IBOUND is negative. 


CARD GROUP 7 (Optional) Format (8E10.0) 


Variable 


Description 


TABLE 


Tables of wall temperature or wall heat flvix versus 

time. IBOUND = -1 these are T versus time tables. 

w 

IBOUND = -2 these are q versus time tables. The 

^w 

input is alternate values of time and T (or q ); three 


w 


*w 


pairs to a card until NPTI pairs are read. (t^,T 


w. 


t„ T , t,, T ). 

2* W 2 3 Wj 


The times are in seconds and the temperatures are in 
®F or °C. The heat fluxes are in watts/cm^ or Btu/ 
ft2-sec. 


CARD GROUP 8 (Optional) Format (El 5.7) 
Variable Description 


TIME Restart time. If ISTART = 2, this card is read, if 

ISTART = 1, do not input it. TIME is the dimensionless 
time at which the program is to be restarted. 


CARD GROUP 9 (Optional) Format (5E15.7) 

Variable Description 

T(I) Temperature profile (I = 1, NN) at the time the run is to 

be restarted (dimensionless). These and the other re- 
start profiles are to be punched from the printout at the 
point where the previous run stopped. 
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CARD GROUP 10 (Optional) Format (5E15.7) 

Variable Description 

U(I) \elocity profile for restart (see Group 9) 

CARD GROUP 11 (Optional) Format (5E15.7) 

Variable Description 

RHO(I) Density profile for restart (see Group 9). 

CARD GROUP 12 Format (80A1) 

Variable Description 

TITLE Run termination card. As many cases as desired can 

be run back-to-back simply by staking them. After the 
last case, a psuedo-title card is punched vsdth the 
characters 

ENDRUN (Cols. 1-6). 

This causes a normal program stop. 

A. 3 SAMPLE CASE 


The infinite parallel plate problem discussed in Section 4 is shown here 
as an example of the input and output formats of the TCI program. A sample 
input is shown below. 

INPUT CARDS 


INFlNlTF Pt ATF PPOFLFM 
20 lonn 1 1 

15.T 273*0 

i.AaoF-^a i*67Ff-a 
t 2.0 0.0 

F\«npi jN 


NFL lUM 
1 1 
1 .PlOF-4 
3. 160F+T 
1 


TW=2 


0.0 

1 .665 


0.05 

1 .0 


1000. n.Q 


1 *C 
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A sample output is given in Fig. A- 3. A printout of the input data is 
shown for checking the user input. A sample flow field is shown at one time 
point to illustrate the format. The meaning of the major variables are the 
same as described in the input guide. 
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Fig. A. 3 — Sample Output of TCI Program 
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DUFORT-FRANKEL NUMERICAL 
SCHEME FOR THERMAL CONVECTION 
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APPENDIX B 

The numerical scheme presented in Section 3 of this report, has first- 
order accuracy truncation error in the time variable. For some fluid mechanics 
applications this order of accuracy is insufficient to produce meaningful results. 
However, for the thermal convection problems being analyzed here, these first- 
order methods should be sufficiently accurate. In order to provide a check on 
the accuracy of the first order method, a Dufort-Frankel (Ref. B-1) algorithm 
was devised and coded. The classical Dufort-Frankel scheme was chosen for 
comparison for the following reasons: 

• Second-order accuracy truncation error in both time 
and distance 

• Simplicity of programming 

• The smallest computer/real time ratio of the five 
schemes tested by Torrance (Ref. B-1) 

The form of the conservation equations and their finite difference repre- 
sentations are shown in Fig. B-1 to B-3. These finite differences incorporate 
the same node-centered spatial mesh used in the original TCI program. 

The TCI explicit technique is utilized as a starter for this multi-time 
step scheme. Solutions are marched out in time from the known initial con- 
ditions. The same logic, schematically accounted for in Appendix A was used 
for the Dufort-Frankel scheme as in TCI except that a CALL START block was 
inserted immediately following the CALL INPUT blbck. 

The Dufort-Frankel program was applied to the sample problem of Larkin 
(Ref. B-2) and the results compared to those of the Lockheed TCI program. 
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Temperature profiles for the sample problenl are compared in Fig. B-4. 

These profiles are practically identical except that the Dufort-Frankel 
temperature gradient at the cool wall is more realistic ® 

rather than positive j . This was a relatively minor change, however, as seen 
in Fig. B-4 The velocity (u) and pressure (P) profiles were also very similar, 

via.. 


TCl Dufort-Frankel 



Period of ^ 

Velocity Waves 10 sec 



sec 


P 

avg 


t = 2000 


1.34 


1.33 


Thus, the TCI and Dufort-Frankel solutions of the sample problem are in very 
close agreement. This is significant, since the two solution techniques use 
entirely different numerical approximations, and both yield similar results. 

The agreement indicates that the Lockheed velocity profiles rather than the 
Larkin (Ref. B-2) and Thuraisamy (Ref. B-3) results, are more realistic, since 
Larkin and Thuraisamy used essentially identical numerical techniques. 
Furthermore, the maximum velocity attained in TCI and D-F solutions seem 
more realistic from a physical viewpoint because Larkin and Thuraisamy' s 
maximum velocities were l/5 the speed of sound whereas those of Lockheed 
were only l/50 of sonic. 

In conclusion, the second-order accurate Dufort-Frankel method has pro- 
duced essentially the same results as the first-order method defined in Section 3. 
The accuracy of the first order method is thus established numerically and has 
been used for all calculations presented in the main text of this report. 
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Continuity Equation 


It = - al- <P“> 


Finite Difference Form 


■ Pi ~ ® ® + U - 0) (pu)“ j 


where 0 = constant, l/Z < 0 < 1. S = At/2 Ax 


Difference Equations 


= Pl+1 "l+l - Oi-1 Vl 


i =2, , . . k-i 


= Pi Uj + P2 U2 


X /«. n n 

‘x<P">k = -'>k"k-'>k-i“k-i 


Resulting Tri-Diagonal System 


A j . T » j . ^ rt 

AiPi., +B.Pi tCjPj^, =D 


X "** i , a««,k 


where 


A c *1+1 

A. = - S u. , 
1 i-l 


e *1+1 

C. e s U., , 
1 1+1 


= 1.0 


D. =p”- S(l-6)Sjj(pu)|‘ 


Aj = 0.0 
Bj = 1 + so 
C, = SO uf * 


= - S 0 u. 


B^ = i-seuf‘ 


S = »•» 


Fig, B-2 - Finite Difference Form of Continuity Equation (Larkin) 





Fig. B-3 - Finite Difference Form of Energy Equation (Dufort-Frankel) 
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